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The intrinsic lattice thermal conductivity of M0S2 is an important aspect in the 
design of MoS2-based nanoelectronic devices. We investigate the lattice dynamics 
properties of M0S2 by first principles calculations. The intrinsic thermal conduc¬ 
tivity of single-layer M0S2 is calculated using the Boltzmann transport equation for 
phonons. The obtained thermal conductivity agrees well with the measurements. 

The contributions of acoustic and optical phonons to the lattice thermal conductiv¬ 
ity are evaluated. The size dependence of thermal conductivity is investigated as 
well. 


I. INTRODUCTION 

Two-dimensional (2D) materials such as graphene and single-layer (SL) transition metal 
dichalcogenides MX 2 (M = Mo, W; X = S, Se, or Te) have drawn considerable interest in 
recent years due to their potential for integration into next-generation electronic and energy 
conversion devices 1 3 . SL M0S2 consists of a Mo atom layer sandwiched between two S atom 
layers, connected by covalent bonds. In the SL form, M0S2 has a direct bandgap of 1.9 eV, 
and switches to an indirect bandgap of 1.6 eV in double-layer M0S2. The bandgap shrinks 
when more layers added, and drops to 1.3 eV in the bulk form as a layered semiconduc¬ 
tor, where the adjacent layers are connected by the van der Waals forcd^. Compared to 
gapless graphene, the extraordinary electronic structure of SL M0S2 may lead to many po¬ 
tential applications such as field-effect transistor^ 6 ( optoelectronic devices 7 ’, and spin-valley 
device J®. 

Since SL M0S2 is considered to play an important role in the next-generation nano¬ 
electronic devices, it is necessary to investigate the thermal properties of M0S2. High- 
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performance electronic devices strongly depend on high thermal conductivity for highly 
efficient heat dissipation, while low thermal conductivity is preferred in thermoelectric appli¬ 
cation. Recently, experimentally measured thermal conductivity k of bulk and SL M 0 S 2 are 
110±20 W/mK 10 and 84±17 W/mK 11 , respectively, since phonons are more sensitive to sur¬ 
face disorder with decreasing layer^ 2 . In addition, theoretically calculated values of n based 
on Boltzmann transport equation (BTE) 13 "^, Molecular-dynamics (MD) simulation^, the 
Klemens model 17 -*®, present great disagreement even by orders of magnitude. Therefore, 
a precise calculation of k of SL M 0 S 2 and clear explanations are needed to clarify such 
disagreements. 

According to the kinetic theory®, accurate prediction of k requires the precise calcu¬ 
lation of the distribution of phonon mean free paths (MFPs). Liu et.al. have measured 
the distribution of phonon MFPs of M 0 S 2 and suggest that approximately 25% of heat is 
carried by phonons with MFP larger than 2 /rm 10 . However, the reported MD simulations 
predicted the longest MFP to be 5.2 nnf 16 *, while some former theoretical calculations based 
on the Klemens’ expression obtained the longest MFP in the scale of 10-200 nm 1 ®, which 
is much lower than the measured value. As a result, the corresponding estimated k is from 
1.35 W/mK to 29.2 W/mK, which is lower than the measured 84±17 W/mK as well. 

In this paper, we obtain phonon properties from lattice dynamics calculations. The 
thermal conductivity of isotopically pure and naturally occurring SL M 0 S 2 are calculated 
using an iterative solution of the BTE for phonons. The calculated thermal conductivity 
agrees well with the measurements. Accurate relaxation time and MFP of M 0 S 2 is obtained. 
The calculated MFP distribution of M 0 S 2 is in consistency to the experimental results. The 
MFP in the small-grain limit is also investigated when the nanostructuring induced phonon 
scattering dominates. The role of boundary scattering in M 0 S 2 nanowires is examined as 
well. 


II. METHODOLOGY 


The in-plane k can be calculated as a sum of contribution of all the phonon modes A, 
which comprises both a phonon branch index p and a wave vector q, 
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where Nk is the number of sampling points, C\ is the heat capacity per mode, V\ a and T\ a 
are the group velocity and relaxation time of mode A along a direction. 

The phonon properties (C\, V \ a , T\ a ) in Eq. Q can be obtained through lattice dy¬ 
namics calculations and an iterative solution of the phonon BTE. The total lattice energy 
can be expanded by a Taylor series with respect to atomic displacements, which includes 
harmonic and (mainly third-order) anharmonic termd- 21 l The harmonic and anharmonic 
interatomic force constants (IFCs) can be obtained from harmonic and anharmonic terms 
respectively. Using the harmonic IFCs, the phonon dispersion relation can be obtained, 
which determines the group velocity v\ a and specific heat C\. The contribution to v\ a and 
C\ from anharmonic terms is usually negligible, and can be thus neglected in the corre¬ 
sponding calculations. However, the third-order term plays an important role in calculating 
the three-phonon scattering rate, which is the inverse of T\ a . 

All the calculations are performed using the Vienna Ab-initio Simulation Package (VASP) 
based on the density functional theory (DFT) 22 . We use the projected augmented wave 
(PAW) method, and the generalized gradient approximation (GGA) in the Perdew-Burke- 
Ernzerhof (PBE) parametrization for the exchange-correlation functional. A cutoff of 500 
eV is used for the plane-wave expansion. A 15x15x1 k-mesh is used for the unit cell during 
structural relaxation. The structures are relaxed until the energy differences are converged 
within 1CU 8 eV, with a Hellman-Feynman force convergence threshold of 10 -4 eV/A. We 
maintain the interlayer vacuum spacing larger than 10 A to eliminate interactions between 
adjacent supercells. 

The harmonic IFCs are obtained by density functional perturbation theory (DFPT) using 
the superccll approach, which calculates the dynamical matrix through the linear response 
of electron density®. A 5x5x1 supercell with 5x5x1 q-mesh is used to calculate the 
dynamical matrix, which is the Fourier transform of the real-space harmonic IFCs < h Q , / g, and 
can be given by 

Aujj(q) = Y y/MjMj ^ a|3 ( Rj ) ex pHq • R A> ( 2 ) 

where a and /3 are the Cartesian indices, and I/ J is the I / J-th atom. The phonon frequen¬ 
cies and eigenvectors can be directly obtained by the solution to the eigenvalue equation 

Y Ay?,/j(q)e/?,j(q) = o; 2 e Q ,/(q). 
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The anharmonic third order IFCs are calculated using a supercell-based, finite-difference 
method^, 

( 4 ) 
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The same 5x5x1 supercell and 5x5x1 q-mesli are used to obtain the anharmonic IFCs. A 
well-converged interaction range of 4.2 A is considered herein, which includes forth-nearest- 
neighbor atoms. 

The lattice thermal conductivity k can be calculated iteratively using the ShengBTE 
code, which is completely parameter-free and based only on the information of the chemical 
structure 24 '" 27 . A discretizationa of the Brillouin zone (BZ) into a T-centered regular grid of 
90x90x1 q points is introduced. 


III. RESULTS AND DISCUSSION 

A. Lattice dynamics properties 

The calculated phonon band structure and phonon density of states (DOS) of SL MoS 2 are 
shown in Fig. [lja) , in which the lattice vibration modes are characterized by three acoustic 
[longitudinal acoustic (LA) and transverse acoustic (TA) branches in the basal plane, and 
flexural acoustic (ZA) branch perpendicular to the basal plane], and six optical branches. 
The LA and TA branches are linear and the ZA branch is quadratic in the vicinity of the 
T point, which is due to the low lattice dimensionality 28 . The calculated bandgap between 
the acoustic and optical branches is about 46 cm -1 , which is in good agreement with other 
previous theoretical result^EI 

According to the group-theoretical analysis, since the SL MoS 2 belongs to the D 3 h point 
group symmetry, the optical lattice-vibration modes at F can be thus decomposed as, 

rgLi = ^'(IR) + 4(R) + UiR + R) + £"(R), (6) 

where IR and R denote infrared- and Raman-active modes respectively. The oscillation 
patterns of the optical modes are shown in Fig. [jjb). Obviously, A!^ and A\ modes are 
out-of-plane vibration modes, while E' and E" are in-plane vibration modes. Table [T] lists 
the calculated frequencies of the four optical phonon modes at the Y point compared to 
experimental values for MoS 2 . The calculated phonon frequencies are in agreement with 
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FIG. 1: (a) Phonon dispersion and DOS of SL M 0 S 2 . (b) Oscillation patterns for phonon modes 
at the T point. 


TABLE I: Phonon frequencies of the four optical modes at the T point compared to experimental 
data (in the units of cm -1 ). 

A '2 A\ E' E" 

This work 458.24 397.63 373.34 276.72 
Experimental 47CP 402.#^ 383. 287® 


the experimental results, and the discrepancy is less than 4%. The LO/TO splitting is very 
small and can be neglected here 17 . 

The Debye temperature do can be calculated with the highest frequency of normal mode 
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FIG. 2: Calculated phonon group velocities of M 0 S 2 in the irreducible BZ. 
vibration (Debye frequency) u m , 

0d = hu m /k B (6) 

where h is the Planck constant, and k B is the Boltzmann constant. The calculated Debye 
temperatures 0 B for M 0 S 2 is 262.3 K, which are in good agreement with previous results, 
i.e. 260-320 K for M 0 S 2 estimated from specific-heat measurement. 33 . The Debye tempera¬ 
ture reflects the magnitude of sound velocity. Higher Debye temperature means increased 
phonon velocities and increased acoustic-phonon frequencies, which suppress phonon-phonon 
scattering by decreasing phonon population^- 34 35 . 

The group velocities v g in the irreducible BZ are shown in Fig [2j which indicate weak 
anisotropy at high frequencies within the entire BZ. The sound velocities in long-wavelength 
limit are about 4,093 and 6,549 m/s for the TA and LA modes respectively, which are 
comparable to 5,400-8,800 m/s in silicent®, and 4,000-8,000 m/s in blue phosphorenc; 3 '. 
The group velocities of optical phonons are less than 3000 m/s, which are much smaller 
than acoustic phonons. 
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FIG. 3 : Lattice thermal conductivity of isotopically pure and naturally occurring SL M0S2, in 
comparison with the results of naturally occurring M0S2 calculated from the single-mode relaxation 
time approximation (SMA). 


B. Intrinsic thermal conductivity 

Fig. [3] presents the lattice thermal conductivity k of isotopically pure and naturally oc¬ 
curring M 0 S 2 using the iterative solution of the BTE. The intrinsic n pure of M 0 S 2 is 101.0 
W/mK, while for naturally occurring isotope concentrations, the K, iso is 87.6 W/ 111 K, which 
agrees well with the measured value of 84±17 W/mK? 11 ! 

Naturally occurring Mo consists of 14.84% 92 Mo, 9.25% 94 Mo, 15.92% 95 Mo, 16.68% 96 Mo, 
9.55% 97 Mo, 24.13% 98 Mo, and 9.63% 100 Mo, while naturally occurring S consists of 95.02% 
32 S, 0.75% 33 S, 4.21% 34 S, and 0.02% 36 S. 

The introduction of isotope disorder can be used to modulate 10% of the thermal conduc¬ 
tivity in isotopically pure M 0 S 2 at temperatures below 262.3 K. The k becomes less sensitive 
to isotopes at high temperatures, since the Umklapp processes (U processes) become frequent 
enough and drive the thermal conductivity at temperatures well above Op 21 38 . 











TABLE II: Contribution of different phonon branches (ZA, TA, LA, and all optical) towards the 
total thermal conductivity in M 0 S 2 , in comparison with graphene and stanene. 

ZA (%) TA (%) LA (%) Optical (%) 


M 0 S 2 

29.1 

30.4 

39.1 

1.4 

graphene 40 

76 

15 

8 

1 

stanene 41 

13.5 

26.9 

57.5 

2.1 


The comparison between k of SL M 0 S 2 with naturally occurring isotope concentrations 
from the single-mode relaxation time approximation (SMA) and the exact solution of the 
BTE as a function of temperature is also shown in Fig. [3j The SMA assumes that individual 
phonon mode is excited independently, which has no memory of the initial phonon distribu¬ 
tion, therefore the SMA is inadequate to describe the momentum-conserving character of the 
Normal processes (N processes) and it works well only if the U processes dominat^ 39 . With 
increasing temperature, the difference between the two approaches declines. Our results 
indicate the N processes play an important role in thermal transport at low temperatures, 
and the U processes become more important with increasing temperature. 

The contributions of different phonon branches to Ki SO are listed in Table [TTJ in compar¬ 
ison with graphene and stanene. It has been reported that the large contribution of ZA 
phonons to the k of graphene is due to a symmetry selection rule in one-atom-thick mate¬ 
rials, which strongly restricts anharmonic phonon-phonon scattering of the ZA phonons 4 ^, 
while the buckled structure in stanene breaks out the out-of-plane symmetry. As for M 0 S 2 , 
ZA and TA phonons contribute almost equally to k, while the LA contribution to k is a bit 
larger than that from the other two phonon modes. Considering the significant difference 
in the contribution of each acoustic phonon mode to the total thermal conductivity in these 
materials, it is worthwhile to perform a detailed investigation of the scattering mechanism 
in M 0 S 2 . 

We extract the frequency-dependent relaxation time of SL M 0 S 2 in Fig. [f](a). At frequen¬ 
cies below 20 cm -1 , the relaxation time of ZA phonons are larger than TA phonons, while 
LA phonons have much short relaxation time. The relaxation time of three acoustic phonon 
modes is comparable with each other at higher frequencies. The frequency-dependent MFP 
of each phonon branch is calculated by l\ = V\T \, as shown in Fig. [4j)b) . The MFP of LA 
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FIG. 4: Frequency dependence of (a) phonon relaxation time and (b) MFP at 300 K. 
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phonons is slightly larger than the other two branches at frequencies above ~25 cm -1 due 
to larger group velocities. Thus the LA phonons contribute the largest part to the total k. 

The calculated longest relaxation time of acoustic phonons (~10 4 ps) is two to three 
orders of magnitude larger than the previous work ( 20-500 psp 19 based on the Klemens’ 
expression^. As a result, the longest MFP in Ref. 17 and 19 (18.1 nm and 221.4 nm, respec¬ 
tively) is much smaller than measured value of >2 pm 10 . In this work, the calculated MFP 
in Fig[i|b) is in consistency to the experimental result. 

C. Size dependence of n 

The MFP distribution plays an important role in describing the behaviour of phonon 
transport within a sample. Generally, when the characteristic length L is smaller than 
the phonon MFP, i.e. L < l\, phonons move ballistically without collisions, and phonon 
transport is in the ballistic regime, in which the thermal conductivity is smaller than the 
prediction by Fourier’s law. When L > l x , phonon transport changes from the ballistic 
regime to the diffusive regime, in which Fourier’s law is valid. The typical procedure to ana¬ 
lyze the MFP distribution is to calculate the cumulative thermal conductivity as a function 

of MFT®33. 

The cumulative thermal conductivity with MFPs below L can be calculated by the fol¬ 
lowing expression, 

«( L ) = ^ Cxv 2 a r Xa . (7) 

A 

The cumulative thermal conductivity with respect to L at 300 K is shown in Fig. [5j( a). The 
accumulation k(L) increases as L increases, until reaching the thermodynamic limit above a 
length Ldiff, which represents the longest mean free path of the heat carriers. The L^ff of 
isotopically pure and naturally occurring M 0 S 2 are 7.5 fin 1 and 7.2 /irn, respectively, which 
are in consistency to the experimental result. 10 . The MFP distribution also suggests that, to 
reduce the lattice thermal conductivity of SL M 0 S 2 , a sample with a characteristic length 
less than ~7 fim is required. 

Furthermore, when the size of sample gets smaller, the nanostructuring-induced phonon 
scattering becomes dominant over the three-phonon scattering, and the small-grain-limit 
reduced k becomes proportional to a constant value lsc ? 16 • We calculate the ratio of the 
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FIG. 5 : (a) Cumulative thermal conductivity as a function of phonon MFP of isotopically pure and 
naturally occurring SL M0S2 at 300 K, in comparison with the results of naturally occurring M0S2 
calculated from the single-mode relaxation time approximation (SMA). (b) Thermal conductivities 
of M0S2 nanowires along [ 100 ] direction as a function of width. 
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thermal conductivity to the thermal conductivity per unit of MFP in the small-grain limit; 
the Isg is found to be 131.5 nm and 114.1 nm for isotopically pure and naturally occurring 
M 0 S 2 at 300 K. This quantity is crucial for the thermal design to modulate the thermal 
conductivity in the small-grain limit, for example nanowires. 

In a nanowire system, phonons with long MFPs will be strongly scattered by the bound¬ 
ary. As a result, the contribution of these phonons to k will be limited. In M 0 S 2 nanowires, 
k decreases with decreasing width, and drops to about half the maximum k in the thermody¬ 
namic limit at widths about 130 nm and 126 nm for isotopically pure and naturally occurring 
M 0 S 2 , as shown in Fig. [5]^b) . Our result indicates that the lattice thermal conductivity of 
M 0 S 2 is sensitive to boundary scattering, and can be further reduced in nanostructures for 
engineering thermal transport in M 0 S 2 . 


IV. CONCLUSION 

We calculate the lattice thermal conductivity k of SL M 0 S 2 using first-principle calcula¬ 
tion and an iterative solution of the BTE for phonons. The introduction of isotopes leads to 
a 10% reduction of k. The intrinsic relaxation time and the distribution of phonon MFPs are 
investigated in detail. The diffusion-limited MFP SL M 0 S 2 is larger than 7 /jrn at 300 K. The 
size dependence of thermal conductivity is investigated as well for the purpose of designing 
nanostructures. Our work provides a fundamental understanding of phonon transport in SL 
M 0 S 2 to predict the thermal performance of MoS 2 -based potential devices. 
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